This is a shortened version of the lab instructions eliminating a few statwide analyses

Part 1 Preliminary Work Setting Up Analysis of Census Block Groups

A few important libraries for spatial analysis are:

install.packages(“maptools”, dependencies = TRUE) install.packages(“spdep”, dependencies = TRUE) install.packages(“leaflet”, dependencies = TRUE) install.packages(“RColorBrewer”, dependencies = TRUE)

library(spdep)
## Loading required package: sp
## Loading required package: Matrix
## Loading required package: spData
## To access larger datasets in this package, install the spDataLarge
## package with: `install.packages('spDataLarge')`
library(maptools)
## Checking rgeos availability: TRUE
library(leaflet)
library(RColorBrewer)

Read the data by setting the directory to point at the right place You will need to download a shp file from Gauchospace and save it in this directory

CA.poly <- readShapePoly(fn = './shp/LPA_Pop_Char_bg.shp')
## Warning: use rgdal::readOGR or sf::st_read

We use the function class to verify the we have a spatial polygons data frame

class(CA.poly)
## [1] "SpatialPolygonsDataFrame"
## attr(,"package")
## [1] "sp"

A SpatialPolygonsDataFrame object brings together the spatial representations of the polygons with data. data: Object of class “data.frame”; attribute table polygons: Object of class “list”; see SpatialPolygons-class plotOrder: Object of class “integer”; see SpatialPolygons-class bbox: Object of class “matrix”; see Spatial-class proj4string: Object of class “CRS”; see CRS-class

The identifying tags of the polygons in the slot are matched with the row names of the data frame to make sure that the correct data rows are associated with the correct spatial object. In this example we have US Census block groups and for each block group we have added the behvior of people and the chractersitics of the households and persons

A Spatial Polygon Data Frame has four components (in the R jargon these called slots)

Component 1: Data contains the variables that are used in the analysis such as number of households with zero cars, vehicle miles of travel. remember these are at the block group level

str(slot(CA.poly, "data"))
summary(slot(CA.poly, "data"))

This shows an output like : ‘data.frame’: 23198 obs. of 105 variables:

Component 2: This is the polygon slot and contains the “shape” information.

The following will create a map of all the polygons representing block groups

plot(CA.poly)

Analysis of the data in these blockgroups

summary(CA.poly)

Component 3: this is the bobox (bounding box of coordinates that is drawn around the boundaries of CA) Component 4: Is the proj4string that contains the projections.

The @ is used to access a specific slot of the spatial data frame

Operations on a column of data is allowed as usual

summary(CA.poly@data$VMT)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##       0   21738   31289   37816   45815 1148168
CA.poly@data$VMTpr = CA.poly@data$VMT/CA.poly@data$n_pr
summary(CA.poly@data$VMTpr)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's 
##    0.00   18.49   22.46   24.06   28.37   90.00      63

First define dummy variables that are based on the type of land use in each block group

CA.poly@data$center = CA.poly@data$LPAgrp == 4
CA.poly@data$suburb = CA.poly@data$LPAgrp == 3
CA.poly@data$exurb = CA.poly@data$LPAgrp == 2
CA.poly@data$rural = CA.poly@data$LPAgrp == 1
CA.poly@data$none = CA.poly@data$LPAgrp == 0

Then estimate a regression model.

VMTprOLS<-lm(VMTpr~suburb+exurb+rural+HHVEH0+HHVEH1+HHVEH2+HHVEH3+HHVEH4+HHVEH5+HHVEH6+HHAGE7, data=CA.poly@data)

VMTprOLS
## 
## Call:
## lm(formula = VMTpr ~ suburb + exurb + rural + HHVEH0 + HHVEH1 + 
##     HHVEH2 + HHVEH3 + HHVEH4 + HHVEH5 + HHVEH6 + HHAGE7, data = CA.poly@data)
## 
## Coefficients:
## (Intercept)   suburbTRUE    exurbTRUE    ruralTRUE       HHVEH0  
##   19.474709     5.599604     5.391621    11.671735    -0.019978  
##      HHVEH1       HHVEH2       HHVEH3       HHVEH4       HHVEH5  
##   -0.009254     0.010847     0.010957     0.046618    -0.147416  
##      HHVEH6       HHAGE7  
##   -0.064597     0.002544
summary(VMTprOLS)
## 
## Call:
## lm(formula = VMTpr ~ suburb + exurb + rural + HHVEH0 + HHVEH1 + 
##     HHVEH2 + HHVEH3 + HHVEH4 + HHVEH5 + HHVEH6 + HHAGE7, data = CA.poly@data)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -30.705  -3.446  -0.754   2.503  63.486 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 19.4747095  0.1044064 186.528  < 2e-16 ***
## suburbTRUE   5.5996042  0.1005532  55.688  < 2e-16 ***
## exurbTRUE    5.3916213  0.1555270  34.667  < 2e-16 ***
## ruralTRUE   11.6717350  0.2011197  58.034  < 2e-16 ***
## HHVEH0      -0.0199781  0.0017925 -11.145  < 2e-16 ***
## HHVEH1      -0.0092536  0.0007676 -12.055  < 2e-16 ***
## HHVEH2       0.0108465  0.0008535  12.708  < 2e-16 ***
## HHVEH3       0.0109570  0.0019188   5.710 1.14e-08 ***
## HHVEH4       0.0466177  0.0040481  11.516  < 2e-16 ***
## HHVEH5      -0.1474161  0.0087433 -16.861  < 2e-16 ***
## HHVEH6      -0.0645967  0.0099456  -6.495 8.47e-11 ***
## HHAGE7       0.0025441  0.0011925   2.133   0.0329 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 5.687 on 23123 degrees of freedom
##   (63 observations deleted due to missingness)
## Multiple R-squared:  0.4257, Adjusted R-squared:  0.4254 
## F-statistic:  1558 on 11 and 23123 DF,  p-value: < 2.2e-16

Test the residuals from this regression model

library(lmtest)
## Loading required package: zoo
## 
## Attaching package: 'zoo'
## The following objects are masked from 'package:base':
## 
##     as.Date, as.Date.numeric
library(sandwich) # I need this for the robust vcov matrix
coeftest(VMTprOLS, vcov = vcovHC(VMTprOLS, type = "const")) # this is the same as in least squares with no White adjsustment
## 
## t test of coefficients:
## 
##                Estimate  Std. Error  t value  Pr(>|t|)    
## (Intercept) 19.47470947  0.10440640 186.5279 < 2.2e-16 ***
## suburbTRUE   5.59960420  0.10055324  55.6880 < 2.2e-16 ***
## exurbTRUE    5.39162132  0.15552699  34.6668 < 2.2e-16 ***
## ruralTRUE   11.67173496  0.20111973  58.0338 < 2.2e-16 ***
## HHVEH0      -0.01997811  0.00179251 -11.1453 < 2.2e-16 ***
## HHVEH1      -0.00925356  0.00076761 -12.0550 < 2.2e-16 ***
## HHVEH2       0.01084651  0.00085354  12.7077 < 2.2e-16 ***
## HHVEH3       0.01095696  0.00191877   5.7104 1.141e-08 ***
## HHVEH4       0.04661772  0.00404807  11.5160 < 2.2e-16 ***
## HHVEH5      -0.14741614  0.00874325 -16.8606 < 2.2e-16 ***
## HHVEH6      -0.06459674  0.00994559  -6.4950 8.472e-11 ***
## HHAGE7       0.00254406  0.00119249   2.1334    0.0329 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
coeftest(VMTprOLS, vcov = vcovHC(VMTprOLS, type = "HC0")) # this is the traditional White's adjustment to the var-Cov of the coefficient estimates
## 
## t test of coefficients:
## 
##                Estimate  Std. Error  t value  Pr(>|t|)    
## (Intercept) 19.47470947  0.10555576 184.4969 < 2.2e-16 ***
## suburbTRUE   5.59960420  0.10105296  55.4126 < 2.2e-16 ***
## exurbTRUE    5.39162132  0.15242269  35.3728 < 2.2e-16 ***
## ruralTRUE   11.67173496  0.20380432  57.2693 < 2.2e-16 ***
## HHVEH0      -0.01997811  0.00209086  -9.5550 < 2.2e-16 ***
## HHVEH1      -0.00925356  0.00084993 -10.8874 < 2.2e-16 ***
## HHVEH2       0.01084651  0.00106403  10.1938 < 2.2e-16 ***
## HHVEH3       0.01095696  0.00211560   5.1791 2.248e-07 ***
## HHVEH4       0.04661772  0.00558563   8.3460 < 2.2e-16 ***
## HHVEH5      -0.14741614  0.01044866 -14.1086 < 2.2e-16 ***
## HHVEH6      -0.06459674  0.00796668  -8.1084 5.383e-16 ***
## HHAGE7       0.00254406  0.00118435   2.1481   0.03172 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# the following are other versions of computing the variance matrix of coefficient estimates
coeftest(VMTprOLS, vcov = vcovHC(VMTprOLS, type = "HC1")) # this is improved White's adjustment 
## 
## t test of coefficients:
## 
##                Estimate  Std. Error  t value  Pr(>|t|)    
## (Intercept) 19.47470947  0.10558315 184.4490 < 2.2e-16 ***
## suburbTRUE   5.59960420  0.10107918  55.3982 < 2.2e-16 ***
## exurbTRUE    5.39162132  0.15246224  35.3637 < 2.2e-16 ***
## ruralTRUE   11.67173496  0.20385720  57.2545 < 2.2e-16 ***
## HHVEH0      -0.01997811  0.00209140  -9.5525 < 2.2e-16 ***
## HHVEH1      -0.00925356  0.00085015 -10.8846 < 2.2e-16 ***
## HHVEH2       0.01084651  0.00106430  10.1912 < 2.2e-16 ***
## HHVEH3       0.01095696  0.00211615   5.1778 2.264e-07 ***
## HHVEH4       0.04661772  0.00558708   8.3438 < 2.2e-16 ***
## HHVEH5      -0.14741614  0.01045138 -14.1050 < 2.2e-16 ***
## HHVEH6      -0.06459674  0.00796875  -8.1063 5.477e-16 ***
## HHAGE7       0.00254406  0.00118465   2.1475   0.03176 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
coeftest(VMTprOLS, vcov = vcovHC(VMTprOLS, type = "HC2")) # this is another improved White's adjustment 
## 
## t test of coefficients:
## 
##                Estimate  Std. Error  t value  Pr(>|t|)    
## (Intercept) 19.47470947  0.10582176 184.0331 < 2.2e-16 ***
## suburbTRUE   5.59960420  0.10120667  55.3284 < 2.2e-16 ***
## exurbTRUE    5.39162132  0.15263767  35.3230 < 2.2e-16 ***
## ruralTRUE   11.67173496  0.20403023  57.2059 < 2.2e-16 ***
## HHVEH0      -0.01997811  0.00209670  -9.5284 < 2.2e-16 ***
## HHVEH1      -0.00925356  0.00085103 -10.8734 < 2.2e-16 ***
## HHVEH2       0.01084651  0.00106576  10.1773 < 2.2e-16 ***
## HHVEH3       0.01095696  0.00212029   5.1677 2.390e-07 ***
## HHVEH4       0.04661772  0.00559303   8.3350 < 2.2e-16 ***
## HHVEH5      -0.14741614  0.01049748 -14.0430 < 2.2e-16 ***
## HHVEH6      -0.06459674  0.00798575  -8.0890 6.309e-16 ***
## HHAGE7       0.00254406  0.00118731   2.1427   0.03215 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
coeftest(VMTprOLS, vcov = vcovHC(VMTprOLS, type = "HC3")) # this is the third improvement
## 
## t test of coefficients:
## 
##                Estimate  Std. Error  t value  Pr(>|t|)    
## (Intercept) 19.47470947  0.10609084 183.5664 < 2.2e-16 ***
## suburbTRUE   5.59960420  0.10136284  55.2432 < 2.2e-16 ***
## exurbTRUE    5.39162132  0.15285505  35.2728 < 2.2e-16 ***
## ruralTRUE   11.67173496  0.20425785  57.1422 < 2.2e-16 ***
## HHVEH0      -0.01997811  0.00210260  -9.5016 < 2.2e-16 ***
## HHVEH1      -0.00925356  0.00085214 -10.8592 < 2.2e-16 ***
## HHVEH2       0.01084651  0.00106751  10.1606 < 2.2e-16 ***
## HHVEH3       0.01095696  0.00212513   5.1559 2.545e-07 ***
## HHVEH4       0.04661772  0.00560063   8.3237 < 2.2e-16 ***
## HHVEH5      -0.14741614  0.01054746 -13.9765 < 2.2e-16 ***
## HHVEH6      -0.06459674  0.00800508  -8.0695 7.401e-16 ***
## HHAGE7       0.00254406  0.00119031   2.1373   0.03258 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
coeftest(VMTprOLS, vcov = vcovHC(VMTprOLS, type = "HC4")) # this is the fourth improvement 
## 
## t test of coefficients:
## 
##                Estimate  Std. Error  t value  Pr(>|t|)    
## (Intercept) 19.47470947  0.10660493 182.6811 < 2.2e-16 ***
## suburbTRUE   5.59960420  0.10164097  55.0920 < 2.2e-16 ***
## exurbTRUE    5.39162132  0.15320295  35.1927 < 2.2e-16 ***
## ruralTRUE   11.67173496  0.20456234  57.0571 < 2.2e-16 ***
## HHVEH0      -0.01997811  0.00211390  -9.4508 < 2.2e-16 ***
## HHVEH1      -0.00925356  0.00085395 -10.8362 < 2.2e-16 ***
## HHVEH2       0.01084651  0.00107056  10.1317 < 2.2e-16 ***
## HHVEH3       0.01095696  0.00213433   5.1337 2.864e-07 ***
## HHVEH4       0.04661772  0.00561357   8.3045 < 2.2e-16 ***
## HHVEH5      -0.14741614  0.01064813 -13.8443 < 2.2e-16 ***
## HHVEH6      -0.06459674  0.00804135  -8.0331 9.957e-16 ***
## HHAGE7       0.00254406  0.00119597   2.1272   0.03341 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
bptest(VMTprOLS, studentize=FALSE )   # the original BP test Page 62 class notes
## 
##  Breusch-Pagan test
## 
## data:  VMTprOLS
## BP = 4989.3, df = 11, p-value < 2.2e-16
bptest(VMTprOLS, studentize=TRUE)     # the studentized version that is more robust to non-normal residuals
## 
##  studentized Breusch-Pagan test
## 
## data:  VMTprOLS
## BP = 1092.9, df = 11, p-value < 2.2e-16

If I run a typical serial correlation test, I may find that there is no correlation between pairs of rows We will check this later

dwtest(VMTprOLS)
## 
##  Durbin-Watson test
## 
## data:  VMTprOLS
## DW = 1.9262, p-value = 9.589e-09
## alternative hypothesis: true autocorrelation is greater than 0

Part 2 Extract Riverside County Data to Perform Analysis of Census Block Groups

We can select a portion of the data. For example, selecting two counties TWOCOUNTY <- CA.poly[CA.poly@data$countyname== c(“Riverside” , “Impreial”), ]

For the class example I just want to work with one county

YCOUNTY <- CA.poly[CA.poly@data$countyname== c("Riverside"), ]

Then I want to know what is in the data slot of this county If you just run YCOUNTY@data you will get all the records - Try just for fun

You can give this instruction to look at the content of the data (not run here to keep output small)

summary(YCOUNTY@data)

Let’s estimate the same regression model we estimated for the entire State of California

VMTprOLSRiver<-lm(VMTpr~suburb+exurb+rural+HHVEH0+HHVEH1+HHVEH2+HHVEH3+HHVEH4+HHVEH5+HHVEH6+HHAGE7, data=YCOUNTY@data)

summary(VMTprOLSRiver)
## 
## Call:
## lm(formula = VMTpr ~ suburb + exurb + rural + HHVEH0 + HHVEH1 + 
##     HHVEH2 + HHVEH3 + HHVEH4 + HHVEH5 + HHVEH6 + HHAGE7, data = YCOUNTY@data)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -17.319  -3.779  -0.688   3.021  49.862 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 21.687987   0.734600  29.524  < 2e-16 ***
## suburbTRUE   5.004393   0.688883   7.265 7.44e-13 ***
## exurbTRUE    2.831428   0.862265   3.284  0.00106 ** 
## ruralTRUE    8.050345   1.153020   6.982 5.25e-12 ***
## HHVEH0      -0.039912   0.009968  -4.004 6.68e-05 ***
## HHVEH1      -0.018997   0.003620  -5.248 1.87e-07 ***
## HHVEH2       0.001059   0.003655   0.290  0.77216    
## HHVEH3       0.020531   0.007350   2.793  0.00532 ** 
## HHVEH4       0.112590   0.015355   7.332 4.61e-13 ***
## HHVEH5      -0.313633   0.038005  -8.253 4.77e-16 ***
## HHVEH6       0.107290   0.056601   1.896  0.05830 .  
## HHAGE7       0.017379   0.003758   4.625 4.23e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 5.685 on 1017 degrees of freedom
##   (1 observation deleted due to missingness)
## Multiple R-squared:  0.4126, Adjusted R-squared:  0.4062 
## F-statistic: 64.93 on 11 and 1017 DF,  p-value: < 2.2e-16

This model is different than what we got using the entire State (see the coefficient for exurbTRUE and compare it to suburbTRUE)

You can also do the BP and DW tests and see what you get but we can do some other more intersting things.

Building Neighborhoods Using Contiguity Rules. The example below uses queen moves (from Chess) to build links among block groups It also creates centroids as centers of “gravity”

In R this function is poly2nb

This builds a neighbours list based on regions with contiguous boundaries. queen=T allows even for a single polygon neighbor to meet the contiguity condition.

list.queenY<-poly2nb(YCOUNTY, queen=T)
coordsY<-coordinates(YCOUNTY)
plot(YCOUNTY)
plot(list.queenY, coordsY, add=T)

summary(list.queenY)
## Neighbour list object:
## Number of regions: 1030 
## Number of nonzero links: 6508 
## Percentage nonzero weights: 0.6134414 
## Average number of links: 6.318447 
## Link number distribution:
## 
##   1   2   3   4   5   6   7   8   9  10  11  12  13  14  15  16  22 
##   1  12  56 102 215 217 186 113  62  28  16   9   3   3   3   2   2 
## 1 least connected region:
## 9993 with 1 link
## 2 most connected regions:
## 20444 20592 with 22 links

Using neighborhoods based on the k-nearest neighbor rule

coords<-coordinates(YCOUNTY)
IDs<-row.names(as(YCOUNTY, "data.frame"))
plot(YCOUNTY)
sids_kn10<-knn2nb(knearneigh(coords, k=10), row.names=IDs)
plot(sids_kn10, coordsY, add=T)

summary (sids_kn10)

House keeping task - Make sure we have valid data in all polygons

One complication is that one of the block groups has NA in VMTpr I change this to zero to avoid missing data and loose a polygon

summary(YCOUNTY@data$VMTpr)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's 
##   11.84   20.66   25.02   25.96   30.48   76.48       1
YCOUNTY@data$VMTpr[is.na(YCOUNTY@data$VMTpr)] <- 0
summary(YCOUNTY@data$VMTpr)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##    0.00   20.63   25.00   25.93   30.48   76.48

Weights creation. This is the matrix W from lecture notes. We create two types of weights Queen and k=10 nearest neighbors weights

Weights without row standardization look like this

$weights[[670]][1] 1 1 1 1 1 1 1

$weights[[671]][1] 1 1 1 1 1 1 1

$weights[[672]][1] 1 1 1 1 1

$weights[[673]][1] 1 1 1 1 1

$weights[[674]][1] 1 1 1

$weights[[675]][1] 1 1 1 1 1 1

$weights[[676]][1] 1 1 1 1 1 1 1 1

Weights with row standardization look like this (I get this out using head(weightsname)

$weights[[679]][1] 0.1428571 0.1428571 0.1428571 0.1428571 0.1428571 0.1428571 0.1428571

$weights[[680]][1] 0.1666667 0.1666667 0.1666667 0.1666667 0.1666667 0.1666667

$weights[[681]][1] 0.1666667 0.1666667 0.1666667 0.1666667 0.1666667 0.1666667

$weights[[682]][1] 0.1666667 0.1666667 0.1666667 0.1666667 0.1666667 0.1666667

$weights[[683]][1] 0.1666667 0.1666667 0.1666667 0.1666667 0.1666667 0.1666667

$weights[[684]][1] 0.25 0.25 0.25 0.25

$weights[[685]][1] 0.2 0.2 0.2 0.2 0.2

$weights[[686]][1] 0.2 0.2 0.2 0.2 0.2

$weights[[687]][1] 0.1111111 0.1111111 0.1111111 0.1111111 0.1111111 0.1111111 0.1111111 0.1111111 0.1111111

Create weights using 0 and 1s for connectivity.

From spdep: nb2listw(neighbours, glist=NULL, style=“codong type of weights”, zero.policy=FALSE)

The function adds a weights list with values given by the coding scheme style chosen. B is the basic binary coding, W is row standardised (sums over all links to n), C is globally standardised (sums over all links to n), U is equal to C divided by the number of neighbours (sums over all links to unity), while S is the variance-stabilizing coding scheme sums over all links to n.

queen_w <- nb2listw(list.queenY, style="B")
summary(queen_w)
## Characteristics of weights list object:
## Neighbour list object:
## Number of regions: 1030 
## Number of nonzero links: 6508 
## Percentage nonzero weights: 0.6134414 
## Average number of links: 6.318447 
## Link number distribution:
## 
##   1   2   3   4   5   6   7   8   9  10  11  12  13  14  15  16  22 
##   1  12  56 102 215 217 186 113  62  28  16   9   3   3   3   2   2 
## 1 least connected region:
## 9993 with 1 link
## 2 most connected regions:
## 20444 20592 with 22 links
## 
## Weights style: B 
## Weights constants summary:
##      n      nn   S0    S1     S2
## B 1030 1060900 6508 13016 184088

Create weights using row standardized weights

queen_ws <- nb2listw(list.queenY)
summary(queen_ws)
## Characteristics of weights list object:
## Neighbour list object:
## Number of regions: 1030 
## Number of nonzero links: 6508 
## Percentage nonzero weights: 0.6134414 
## Average number of links: 6.318447 
## Link number distribution:
## 
##   1   2   3   4   5   6   7   8   9  10  11  12  13  14  15  16  22 
##   1  12  56 102 215 217 186 113  62  28  16   9   3   3   3   2   2 
## 1 least connected region:
## 9993 with 1 link
## 2 most connected regions:
## 20444 20592 with 22 links
## 
## Weights style: W 
## Weights constants summary:
##      n      nn   S0       S1     S2
## W 1030 1060900 1030 345.3927 4291.3

The same weight creation but using k nearest neighbor

sids_kn10_w<- nb2listw(sids_kn10)
summary(sids_kn10_w)
## Characteristics of weights list object:
## Neighbour list object:
## Number of regions: 1030 
## Number of nonzero links: 10300 
## Percentage nonzero weights: 0.9708738 
## Average number of links: 10 
## Non-symmetric neighbours list
## Link number distribution:
## 
##   10 
## 1030 
## 1030 least connected regions:
## 9 101 107 134 174 318 364 420 454 493 501 520 525 533 546 563 575 628 642 670 695 699 724 742 821 822 886 1130 1198 1233 1239 1279 1290 1295 1312 1345 1400 1405 1411 1441 1454 1506 1527 1553 1556 1609 1650 1656 1691 1697 1703 1714 1744 1746 1790 1813 1844 1884 1900 1905 1981 1995 2033 2079 2104 2249 2253 2270 2289 2298 2326 2352 3964 3965 3966 3967 3968 3969 3970 3971 3972 3973 3974 3975 3976 3977 3978 3979 3980 3981 3982 3983 3984 3985 3986 3987 3988 3989 4001 4002 4003 4004 4005 4111 4112 4113 4114 4115 4116 4117 4118 4119 4120 4121 4122 4123 4124 4125 4126 4127 4128 4129 4130 4131 4132 4133 4134 4135 4136 4137 4138 4139 4140 4141 4142 4143 4144 4145 4146 4147 4148 4149 4150 4151 4702 4703 4704 4705 4706 4707 4708 4709 4710 4711 4712 4713 4714 4715 4716 4717 4718 4719 4720 4721 4722 4723 4724 4725 4726 4727 4728 4729 4730 4731 4732 4733 4734 4735 4736 4737 4738 4739 4740 4741 4742 4743 4744 4745 4746 4747 4748 4749 4750 4751 4752 4753 4754 4755 4756 4757 4758 4759 4760 4761 4762 4763 4764 4765 4766 4767 4768 4769 4770 4771 4772 4773 4774 4775 4776 4777 4778 4779 4780 5356 5412 5419 5443 5509 5516 5563 5609 5735 5739 5872 5911 5947 5992 6081 6082 6103 6116 6128 6142 6153 6194 6235 6253 6264 6268 6270 6278 6296 6344 6369 6379 6382 6428 6437 6457 6466 6467 6490 6492 6509 6513 6514 6546 6549 6587 6599 6603 6655 6702 6736 6791 6851 6864 6886 7024 7036 7176 7222 7260 7277 7295 7318 7400 7404 7510 7522 7801 8173 8174 8175 8176 8177 8178 8179 8184 8185 8186 8187 8188 8189 8190 8197 8198 8199 8200 8201 8202 8211 8212 8213 8214 8436 8437 8438 8439 8440 8441 8442 8443 8449 8450 8451 8452 8453 8454 8455 8456 8457 8458 8459 8462 8463 8466 8467 8468 8469 8470 8657 8660 8730 8747 8775 8800 8802 8896 8922 8930 8947 9032 9033 9035 9050 9099 9194 9199 9221 9224 9236 9257 9593 9607 9608 9609 9610 9623 9672 9673 9674 9679 9680 9681 9691 9692 9697 9698 9699 9700 9701 9702 9703 9704 9705 9706 9818 9866 9899 9932 9934 9993 10006 10026 10029 10033 10036 10049 10069 10100 10170 10180 10193 10198 10289 10345 10394 10490 10491 10492 10493 10495 10525 10526 10531 10532 10533 10534 10536 10537 10538 10653 10661 10663 10672 10690 10704 10714 10754 10775 10794 10838 10863 10898 10900 10906 10907 10908 10910 10924 10928 11013 11044 11046 11071 11098 11112 11146 11158 11170 11172 11217 11232 11236 11242 11290 11292 11297 11303 11329 11332 11336 11344 11365 11371 11379 11388 11393 11429 11445 11520 11522 11527 11531 11552 11553 11554 11562 11629 11638 11701 11717 11719 11722 11769 11792 11803 11804 11825 11839 11894 11917 11922 11966 11985 11992 12030 12045 12069 12073 12122 12133 12145 12149 12251 12267 12288 12302 12306 12312 12344 12365 12406 12534 12542 12572 12608 12642 12651 12663 12681 12690 12701 12706 12708 12766 12775 12780 12915 12916 12928 12947 12953 13004 13034 13037 13056 13091 13096 13120 13128 13140 13148 13174 13179 13277 13291 13315 13364 13371 13390 13394 13426 13428 13448 13511 13548 13565 13730 13735 15531 15532 15533 15534 15535 15536 15537 15538 15539 15547 15548 15549 15550 15551 15565 15566 15567 15568 15569 15570 15571 15572 15573 15574 15575 15576 15577 15589 15590 15591 15592 15593 15594 15595 15596 15617 15618 15635 15636 15637 15638 15639 15640 15688 15689 15690 15691 15692 15707 15708 15709 15710 15711 15712 15713 15714 15715 15761 15762 15763 15764 15765 15766 15767 15768 15769 15770 15771 15772 15773 15774 15775 15796 15797 15798 15799 15800 15801 15802 15803 15804 15805 15843 15844 15845 15846 15847 15848 15849 15850 15860 15861 15862 15863 15864 15886 15887 16511 16512 16534 16535 16536 16537 16538 16539 16540 16541 16542 16543 16544 16545 16546 16547 16548 16549 16550 16551 16552 16553 16554 16555 16556 16557 16558 16559 16560 16561 16562 16563 16564 16565 16566 16567 16568 16569 16570 16571 16572 16573 16574 16575 16576 16577 16578 16579 16580 16581 16582 16583 16584 16585 16586 16587 16588 16589 16590 16591 16592 16593 16594 16595 16596 16597 16598 16599 16600 16601 16602 16603 16604 16605 16606 16607 16608 16609 16610 16611 16612 16613 16614 16615 16616 16617 16618 16619 16620 16621 16622 16623 16624 16625 16626 16627 16628 17383 17411 17424 17429 17518 17521 17525 17535 17572 17587 17600 17621 17625 17633 17704 17705 17708 17720 17721 17728 17731 17732 17787 17790 17812 17817 17830 17848 17913 17916 17932 17943 18005 18007 18035 18045 18059 18065 18097 18114 18151 18170 18249 18260 18296 18307 18339 18371 18386 18413 18459 18484 18490 18502 18523 18532 18543 18556 18558 18618 18651 18674 18677 18697 18719 18776 18793 18858 18876 18909 18931 18944 18970 19044 19055 19056 19085 19099 19134 19142 19147 19172 19203 19208 19226 19247 19339 19353 19356 19372 19374 19392 19409 19411 19439 19447 19481 19520 19536 19581 19672 19698 19783 19786 19787 19791 19837 19958 20132 20357 20398 20415 20416 20417 20418 20419 20420 20421 20444 20516 20517 20518 20519 20520 20521 20544 20545 20546 20556 20566 20567 20568 20588 20592 20593 20594 20595 20656 20679 20680 20681 20686 20687 20688 20689 20690 20691 20726 20727 20728 20729 20730 20731 20732 20733 20734 20735 20739 20740 20741 20742 20743 20774 20775 20776 20777 20778 20779 20782 20794 20795 20796 20797 21011 21078 21079 21093 21176 21191 21277 21282 21294 21321 21326 21332 21359 21444 21447 21449 21463 21477 21509 21525 21530 21588 21607 21695 21727 21922 21926 21927 21928 21930 21931 21932 21935 21943 21944 21955 21970 21971 21972 21983 21984 22008 22022 22023 22024 22025 22026 22027 22028 22037 22038 22039 22040 22041 22042 22064 22083 22094 22177 22242 22300 22311 22323 22360 22428 22443 22450 22456 22462 22496 22505 22532 22564 22682 22683 22690 22743 22778 22815 22816 22817 22818 22825 22871 22872 22888 22889 22913 22923 22924 22925 22953 22988 23001 23011 23019 23107 23113 23116 23120 23164 23165 with 10 links
## 1030 most connected regions:
## 9 101 107 134 174 318 364 420 454 493 501 520 525 533 546 563 575 628 642 670 695 699 724 742 821 822 886 1130 1198 1233 1239 1279 1290 1295 1312 1345 1400 1405 1411 1441 1454 1506 1527 1553 1556 1609 1650 1656 1691 1697 1703 1714 1744 1746 1790 1813 1844 1884 1900 1905 1981 1995 2033 2079 2104 2249 2253 2270 2289 2298 2326 2352 3964 3965 3966 3967 3968 3969 3970 3971 3972 3973 3974 3975 3976 3977 3978 3979 3980 3981 3982 3983 3984 3985 3986 3987 3988 3989 4001 4002 4003 4004 4005 4111 4112 4113 4114 4115 4116 4117 4118 4119 4120 4121 4122 4123 4124 4125 4126 4127 4128 4129 4130 4131 4132 4133 4134 4135 4136 4137 4138 4139 4140 4141 4142 4143 4144 4145 4146 4147 4148 4149 4150 4151 4702 4703 4704 4705 4706 4707 4708 4709 4710 4711 4712 4713 4714 4715 4716 4717 4718 4719 4720 4721 4722 4723 4724 4725 4726 4727 4728 4729 4730 4731 4732 4733 4734 4735 4736 4737 4738 4739 4740 4741 4742 4743 4744 4745 4746 4747 4748 4749 4750 4751 4752 4753 4754 4755 4756 4757 4758 4759 4760 4761 4762 4763 4764 4765 4766 4767 4768 4769 4770 4771 4772 4773 4774 4775 4776 4777 4778 4779 4780 5356 5412 5419 5443 5509 5516 5563 5609 5735 5739 5872 5911 5947 5992 6081 6082 6103 6116 6128 6142 6153 6194 6235 6253 6264 6268 6270 6278 6296 6344 6369 6379 6382 6428 6437 6457 6466 6467 6490 6492 6509 6513 6514 6546 6549 6587 6599 6603 6655 6702 6736 6791 6851 6864 6886 7024 7036 7176 7222 7260 7277 7295 7318 7400 7404 7510 7522 7801 8173 8174 8175 8176 8177 8178 8179 8184 8185 8186 8187 8188 8189 8190 8197 8198 8199 8200 8201 8202 8211 8212 8213 8214 8436 8437 8438 8439 8440 8441 8442 8443 8449 8450 8451 8452 8453 8454 8455 8456 8457 8458 8459 8462 8463 8466 8467 8468 8469 8470 8657 8660 8730 8747 8775 8800 8802 8896 8922 8930 8947 9032 9033 9035 9050 9099 9194 9199 9221 9224 9236 9257 9593 9607 9608 9609 9610 9623 9672 9673 9674 9679 9680 9681 9691 9692 9697 9698 9699 9700 9701 9702 9703 9704 9705 9706 9818 9866 9899 9932 9934 9993 10006 10026 10029 10033 10036 10049 10069 10100 10170 10180 10193 10198 10289 10345 10394 10490 10491 10492 10493 10495 10525 10526 10531 10532 10533 10534 10536 10537 10538 10653 10661 10663 10672 10690 10704 10714 10754 10775 10794 10838 10863 10898 10900 10906 10907 10908 10910 10924 10928 11013 11044 11046 11071 11098 11112 11146 11158 11170 11172 11217 11232 11236 11242 11290 11292 11297 11303 11329 11332 11336 11344 11365 11371 11379 11388 11393 11429 11445 11520 11522 11527 11531 11552 11553 11554 11562 11629 11638 11701 11717 11719 11722 11769 11792 11803 11804 11825 11839 11894 11917 11922 11966 11985 11992 12030 12045 12069 12073 12122 12133 12145 12149 12251 12267 12288 12302 12306 12312 12344 12365 12406 12534 12542 12572 12608 12642 12651 12663 12681 12690 12701 12706 12708 12766 12775 12780 12915 12916 12928 12947 12953 13004 13034 13037 13056 13091 13096 13120 13128 13140 13148 13174 13179 13277 13291 13315 13364 13371 13390 13394 13426 13428 13448 13511 13548 13565 13730 13735 15531 15532 15533 15534 15535 15536 15537 15538 15539 15547 15548 15549 15550 15551 15565 15566 15567 15568 15569 15570 15571 15572 15573 15574 15575 15576 15577 15589 15590 15591 15592 15593 15594 15595 15596 15617 15618 15635 15636 15637 15638 15639 15640 15688 15689 15690 15691 15692 15707 15708 15709 15710 15711 15712 15713 15714 15715 15761 15762 15763 15764 15765 15766 15767 15768 15769 15770 15771 15772 15773 15774 15775 15796 15797 15798 15799 15800 15801 15802 15803 15804 15805 15843 15844 15845 15846 15847 15848 15849 15850 15860 15861 15862 15863 15864 15886 15887 16511 16512 16534 16535 16536 16537 16538 16539 16540 16541 16542 16543 16544 16545 16546 16547 16548 16549 16550 16551 16552 16553 16554 16555 16556 16557 16558 16559 16560 16561 16562 16563 16564 16565 16566 16567 16568 16569 16570 16571 16572 16573 16574 16575 16576 16577 16578 16579 16580 16581 16582 16583 16584 16585 16586 16587 16588 16589 16590 16591 16592 16593 16594 16595 16596 16597 16598 16599 16600 16601 16602 16603 16604 16605 16606 16607 16608 16609 16610 16611 16612 16613 16614 16615 16616 16617 16618 16619 16620 16621 16622 16623 16624 16625 16626 16627 16628 17383 17411 17424 17429 17518 17521 17525 17535 17572 17587 17600 17621 17625 17633 17704 17705 17708 17720 17721 17728 17731 17732 17787 17790 17812 17817 17830 17848 17913 17916 17932 17943 18005 18007 18035 18045 18059 18065 18097 18114 18151 18170 18249 18260 18296 18307 18339 18371 18386 18413 18459 18484 18490 18502 18523 18532 18543 18556 18558 18618 18651 18674 18677 18697 18719 18776 18793 18858 18876 18909 18931 18944 18970 19044 19055 19056 19085 19099 19134 19142 19147 19172 19203 19208 19226 19247 19339 19353 19356 19372 19374 19392 19409 19411 19439 19447 19481 19520 19536 19581 19672 19698 19783 19786 19787 19791 19837 19958 20132 20357 20398 20415 20416 20417 20418 20419 20420 20421 20444 20516 20517 20518 20519 20520 20521 20544 20545 20546 20556 20566 20567 20568 20588 20592 20593 20594 20595 20656 20679 20680 20681 20686 20687 20688 20689 20690 20691 20726 20727 20728 20729 20730 20731 20732 20733 20734 20735 20739 20740 20741 20742 20743 20774 20775 20776 20777 20778 20779 20782 20794 20795 20796 20797 21011 21078 21079 21093 21176 21191 21277 21282 21294 21321 21326 21332 21359 21444 21447 21449 21463 21477 21509 21525 21530 21588 21607 21695 21727 21922 21926 21927 21928 21930 21931 21932 21935 21943 21944 21955 21970 21971 21972 21983 21984 22008 22022 22023 22024 22025 22026 22027 22028 22037 22038 22039 22040 22041 22042 22064 22083 22094 22177 22242 22300 22311 22323 22360 22428 22443 22450 22456 22462 22496 22505 22532 22564 22682 22683 22690 22743 22778 22815 22816 22817 22818 22825 22871 22872 22888 22889 22913 22923 22924 22925 22953 22988 23001 23011 23019 23107 23113 23116 23120 23164 23165 with 10 links
## 
## Weights style: W 
## Weights constants summary:
##      n      nn   S0     S1      S2
## W 1030 1060900 1030 180.98 4232.62

I you use style=“C” gives equal weights to all connections and produces the following in terms of weights

$weights[[671]][1] 0.1582667 0.1582667 0.1582667 0.1582667 0.1582667 0.1582667 0.1582667

$weights[[672]][1] 0.1582667 0.1582667 0.1582667 0.1582667 0.1582667

$weights[[673]][1] 0.1582667 0.1582667 0.1582667 0.1582667 0.1582667

$weights[[674]][1] 0.1582667 0.1582667 0.1582667

Autocorrelations at different lags using the defaults in sp.correlogram.

The question we ask is: Do we see spatial correlation in the vehicle miles per person produce by each Censusblock groups? How far is this correlation extending? What difference does it make to use weights using Queen continguity vs k-nearest neighbor continguity?

The following are the plots from slides 57 and 58 of class notes

mor10q <- sp.correlogram(list.queenY, var=YCOUNTY@data$VMTpr, order=10, method="I")
plot(mor10q, main = "Moran's I with Queen Contiguity and Row Standardization")

mor10k <- sp.correlogram(sids_kn10, var=YCOUNTY@data$VMTpr, order=10, method="I", zero.policy=TRUE)  # need zero policy because some polygons are not connected
plot(mor10k, main = "Moran I with knn10 Contiguity and Row Standardization")

ger10q <- sp.correlogram(list.queenY, var=YCOUNTY@data$VMTpr, order=10, method="C")
plot(ger10q, main = "Geary's C with Queen Contiguity and Row Standardization")

ger10k <- sp.correlogram(sids_kn10, var=YCOUNTY@data$VMTpr, order=10, method="C", zero.policy=TRUE)  # need zero policy because some polygons are not connected
plot(ger10k, main = "Geary's C with knn10 Contiguity and Row Standardization")

We can also create statstical tests for spatial correlation based on the z-scores

The first uses a weights matrix that is derived from Queen continguity using the default spatial weights row standardization and computing the variance with analytical randomization

moran.test(YCOUNTY@data$VMTpr, listw=nb2listw(list.queenY))
## 
##  Moran I test under randomisation
## 
## data:  YCOUNTY@data$VMTpr  
## weights: nb2listw(list.queenY)  
## 
## Moran I statistic standard deviate = 14.137, p-value < 2.2e-16
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##      0.2529411655     -0.0009718173      0.0003225841

In the Table above

standarddeviate=((0.2529411655-(-0.0009718173))/sqrt(0.0003225841))
standarddeviate
## [1] 14.1372

This indicates that we have a strong spatial correlation

We repeat the same but this time we assume normality

moran.test(YCOUNTY@data$VMTpr, listw=nb2listw(list.queenY), randomisation=FALSE)
## 
##  Moran I test under normality
## 
## data:  YCOUNTY@data$VMTpr  
## weights: nb2listw(list.queenY)  
## 
## Moran I statistic standard deviate = 14.117, p-value < 2.2e-16
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##      0.2529411655     -0.0009718173      0.0003235222

The standard deviate is not very different from the previous result.
This usually happens when we have many units (the polygons)

The next question is: are there observations that have very high spatial correlations?

spdep has a function called moran.plot. This produces an object with some useful quantities For example it runs a regression of x (the variable we analyze) on wx (the weighted values of the neighbors of x) It also plots each x and wx pairs and the mean of x and wx

Below we use Queen continuity and style=“C” This weights up observations with many neighbors because all connections to all polygons take the same value. Polygons with 6 connections will be influneced by 2 more neighbors than polygons with 4 connections (all weighted with the same value (0.15… we saw before))

msp <- moran.plot(YCOUNTY@data$VMTpr, listw=nb2listw(list.queenY, style="C"), quiet=TRUE)
title("Moran scatterplot Riverside")

The following code is from Bivand/Pebesma/Gomez-Rubio and extracts from object msp the influence of each pair x and wx in the regression of wx on x. This is one way to identify which polygons influences more the global Moran’s index

It also plots a map

infl <- apply(msp$is.inf, 1, any)
x <- YCOUNTY@data$VMTpr
lhx <- cut(x, breaks=c(min(x), mean(x), max(x)), labels=c("L", "H"), include.lowest=TRUE)
wx <- lag(nb2listw(list.queenY, style="C"), YCOUNTY@data$VMTpr)
lhwx <- cut(wx, breaks=c(min(wx), mean(wx), max(wx)), labels=c("L", "H"), include.lowest=TRUE)
lhlh <- interaction(lhx, lhwx, infl, drop=TRUE)
cols <- rep(1, length(lhlh))
cols[lhlh == "None"] <- 1
cols[lhlh == "H.L.TRUE"] <- 2
cols[lhlh == "L.H.TRUE"] <- 3
cols[lhlh == "H.H.TRUE"] <- 4
plot(YCOUNTY, col=brewer.pal(4, "Accent")[cols])
# RSB quietening greys
legend("topright", legend=c("None", "HL", "LH", "HH"), fill=brewer.pal(4, "Accent"), bty="n", cex=0.8, y.intersp=0.8)
title("Block groups with influence")

Not very nice and cannot even tell where the polygons (blockgroups). Leaflet will do the job.

First I define the colors and labels for the categories of influence and then build a map

colsF <- factor(cols, labels=c("None", "HighLow", "LowHigh", "HighHigh"))
LHpal <- colorFactor(topo.colors(4), colsF)
popup <- paste0("<strong> BLOCKGROUP </strong>", IDs)
leaflet(YCOUNTY) %>%
  addPolygons(stroke = FALSE, smoothFactor = 0.2, fillOpacity = 0.7, popup=popup,
              color = ~LHpal(colsF)) %>% addTiles() %>% addLegend("topright", pal = LHpal, values = ~colsF,
                                                                   title = "Influence",
                                                                   opacity = 1)

The definiton weights has a substantial impact on the influence of neighbors and ultimately on Moran’s I From spdep: B is the basic binary coding, W is row standardised (sums over all links to n), C is globally standardised (sums over all links to n), U is equal to C divided by the number of neighbours (sums over all links to unity), while S is the variance-stabilizing coding scheme sums over all links to n.

mspdefault <- moran.plot(YCOUNTY@data$VMTpr, listw=nb2listw(list.queenY), quiet=TRUE)
title("Moran scatterplot Riverside Default Style Default")

mspW <- moran.plot(YCOUNTY@data$VMTpr, listw=nb2listw(list.queenY, style="W"), quiet=TRUE)
title("Moran scatterplot Riverside  Style W")

mspC <- moran.plot(YCOUNTY@data$VMTpr, listw=nb2listw(list.queenY, style="C"), quiet=TRUE)
title("Moran scatterplot Riverside  Style C")

mspB <- moran.plot(YCOUNTY@data$VMTpr, listw=nb2listw(list.queenY, style="B"), quiet=TRUE)
title("Moran scatterplot Riverside  Style B")

mspU <- moran.plot(YCOUNTY@data$VMTpr, listw=nb2listw(list.queenY, style="U"), quiet=TRUE)
title("Moran scatterplot Riverside  Style U")

mspS <- moran.plot(YCOUNTY@data$VMTpr, listw=nb2listw(list.queenY, style="S"), quiet=TRUE)
title("Moran scatterplot Riverside  Style S")

The equation for local Moran’s I (one of the Anselin LISA indicators) - see Gauchospace paper is below. In this book Bivan et al show the denominator division by n not n-1. Numerically does not matter. Division by n is the defaul in spdep

\[ I_i = \frac{(x_i-\bar{x})}{{∑_{k=1}^{n}(x_k-\bar{x})^2}/(n-1)}{∑_{j=1}^{n}w_{ij}(x_j-\bar{x})} \]

We get as many of these indicators as the number of units (blockgroups in our example)

localM1 <- as.data.frame(localmoran(YCOUNTY@data$VMTpr, listw=nb2listw(list.queenY, style="C")))
summary(localM1)
##        Ii                E.Ii                Var.Ii       
##  Min.   :-2.26283   Min.   :-0.0033837   Min.   :0.02493  
##  1st Qu.:-0.02368   1st Qu.:-0.0010766   1st Qu.:0.12415  
##  Median : 0.09788   Median :-0.0009228   Median :0.14884  
##  Mean   : 0.24063   Mean   :-0.0009718   Mean   :0.15658  
##  3rd Qu.: 0.42794   3rd Qu.:-0.0007690   3rd Qu.:0.17348  
##  Max.   : 4.69817   Max.   :-0.0001538   Max.   :0.53725  
##       Z.Ii            Pr(z > 0)     
##  Min.   :-5.43029   Min.   :0.0000  
##  1st Qu.:-0.05803   1st Qu.:0.1280  
##  Median : 0.24817   Median :0.4020  
##  Mean   : 0.61887   Mean   :0.3710  
##  3rd Qu.: 1.13600   3rd Qu.:0.5231  
##  Max.   :11.78184   Max.   :1.0000

Store the Local Moran I and its zscores in the database

YCOUNTY@data$localM1 <- localM1[,1]
summary(YCOUNTY@data$localM1)
##     Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
## -2.26283 -0.02368  0.09788  0.24063  0.42794  4.69817
YCOUNTY@data$zscoreM1 <- localM1[,4]
summary(YCOUNTY@data$zscoreM1)
##     Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
## -5.43029 -0.05803  0.24817  0.61887  1.13600 11.78184

Map the observed x

qpalreds <- colorNumeric(
  palette = "Reds",
  domain = CA.poly@data$VMTpr)

leaflet(YCOUNTY) %>%
  addPolygons(stroke = FALSE, fillOpacity = .8, smoothFactor = 0.2, 
              color = ~qpalreds(VMTpr)) %>% addTiles() %>% addLegend("topright", pal = qpalreds, values = ~VMTpr,
                                                                     title = "Observed Miles per Person",
                                                                     labFormat = labelFormat(prefix = "Observed Miles per person day "),
                                                                     opacity = 0.8)
qpalM <- colorNumeric(
  palette = "Blues",
  domain = CA.poly@data$localM1)

leaflet(YCOUNTY) %>%
  addPolygons(stroke = FALSE, fillOpacity = .8, smoothFactor = 0.2, 
              color = ~qpalM(localM1)) %>% addTiles() %>% addLegend("topright", pal = qpalM, values = ~localM1,
                                                                     title = "Local Moran I",
                                                                     labFormat = labelFormat(prefix = "Local Moran I "),
                                                                     opacity = 1.0)
qpalZ <- colorNumeric(
  palette = "Greens",
  domain = CA.poly@data$zscoreM1)

leaflet(YCOUNTY) %>%
  addPolygons(stroke = FALSE, fillOpacity = .8, smoothFactor = 0.2, 
              color = ~qpalZ(zscoreM1)) %>% addTiles() %>% addLegend("topright", pal = qpalZ, values = ~zscoreM1,
                                                                       title = "Local Moran I Z score",
                                                                       labFormat = labelFormat(prefix = "Local Moran I Z score "),
                                                                       opacity = 0.9)